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User Guide for Compressible Flow Toolbox 

Version 2.1 for Use With MATLAB h Version 7 
Kevin J. Melcher 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 


Abstract 

This report provides a user guide for the Compressible Flow Toolbox, a collection of algorithms 
that solve almost 300 linear and nonlinear classical compressible flow relations. The algorithms, 
implemented in the popular MATLAB K programming language, are useful for analysis of one- 
dimensional steady flow with constant entropy, friction, heat transfer, or shock discontinuities. 
The solutions do not include any gas dissociative effects. The toolbox also contains functions for 
comparing and validating the equation-solving algorithms against solutions previously published 
in the open literature. The classical equations solved by the Compressible Flow Toolbox are: 
isentropic-flow equations, Fanno flow equations (pertaining to flow of an ideal gas in a pipe with 
friction), Rayleigh flow equations (pertaining to frictionless flow of an ideal gas, with heat 
transfer, in a pipe of constant cross section.), normal-shock equations, oblique-shock equations, 
and Prandtl-Meyer expansion equations. At the time this report was published, the Compressible 
Flow Toolbox was available without cost from the NASA Software Repository. 

1. Introduction 

Description 

This paper provides a User Guide for the Compressible Flow Toolbox, a collection of algorithms 
that solve almost 300 linear and nonlinear classical compressible flow relations. The algorithms, 
implemented in the popular MATLAB R programming language, are useful for analysis of one- 
dimensional steady flow with constant entropy, friction, heat transfer, or shock discontinuities. 
The solutions do not include any gas dissociative effects. The toolbox also contains functions for 
comparing and validating the equation-solving algorithms against solutions previously published 
in the open literature. The classical equations solved by the Compressible Flow Toolbox are: 

• The isentropic-flow equations, 

• The Fanno flow equations (pertaining to flow of an ideal gas in a pipe with friction), 

• The Rayleigh flow equations (pertaining to frictionless flow of an ideal gas, with heat 
transfer, in a pipe of constant cross section.) 

• The normal-shock equations, 

• The oblique-shock equations, and 

• The Prandtl-Meyer expansion equations. 

The user should note that the scope of this guide is limited to documenting the individual 
functions and providing instruction in using them to solve simple compressible flow examples. 
Functions in the toolbox can be used together to solve more complex compressible flow 
problems — that is why they were created. However, instructing the user in the broader context of 
compressible flow is not the intended purpose of this guide. 
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Background 

Algorithms included in the Compressible Flow Toolbox were originally developed to support 
controls and dynamics research under the NASA’s High Speed Research Program. They were 
inspired by NACA Report 1135 “Equations Tables, and Charts for Compressible Flow” (ref. 1) 
which the author studied extensively as part of that research. Early implementations were first 
published as part of the author’s Masters Thesis in 1996. They were subsequently made publicly 
available via a M ATLAB ’ third party software web site hosted by the Mathworks, Inc. After 
several years, the toolbox was removed from the web site for a variety of reasons, including the 
need to upgrade the algorithms for compatibility with newer versions of MATLAB ” . Finally, to 
appease a number of recent requests for the software, the toolbox has been updated, expanded, 
and made available to the general public via the NASA Software Repository. 

All of the numerical and graphical results shown in this report were generated using functions 
included in the Compressible Flow Toolbox version 2. land running M ATLAB ’ version 7.04 on 
an MS Windows XP, 2.2 GHz Intel Pentium 4 processor-based personal computer. Results may 
vary slightly based on the precision of the floating point processor used to perform the 
calculations. 

Organization 

This User’s Guide is organized in five sections. Introduction, Nomenclature, Quick Reference 
Guide, Function Reference Guide, and References. Section 1. Introduction provides a general 
description of the User Guide along with historical information on the origin of the toolbox and 
availability of the software. Section 2. Nomenclature describes the symbols and special 
formatting conventions used throughout the text. Section 3. Quick Reference Guide provides a 
comprehensive list of the functions contained in the toolbox and provides a brief description of 
each function listed. Section 4, Function Reference Guide provides a detailed description of 
each function in the toolbox including its purpose, syntax, a discussion of how the algorithm 
works, and examples demonstrating its use. Finally, Section 5. References contains a list of 
references used in developing and documenting the toolbox. 

Availability 

At the time this report was published, the Compressible Flow Toolbox was available to the 
general public without cost through the NASA Software Repository. 

https://technology.grc.nasa.gov/software/ 
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2. Nomenclature 


Formats and Convensions 

Monospace MATLAB® commands, functions names and screen output are 

displayed in this font. For example: rayl ei gh. 


Italics Book titles and names of book sections, mathematical symbols 

and notation, and the introduction of new terms. For example: 
Introduction. 


Bold Initial Caps Key names, menu names, and items that are selected from 
menus. For example: the File menu. 


Symbols 

This document uses the following symbols and notations: 


Roman Symbols 

A 

D h 

I 

M 

P 

P, 

T 

T, 

V 


Description 

Cross-sectional area of stream tube or channel 

Hydraulic diameter of the flow cross-sectional area 

Impulse function 

Mach number, Via 

Static Pressure 

Total Pressure 

Static Temperature 

Total Temperature 

Velocity 


7 

q 


Average friction factor 
Dynamic pressure, p V 2 12 


Greek Symbols Description 


P 

y 

8 

0 

F 

v 

P 

P 


Ratio of specific heats of the working fluid ( default = 1.4) 
Turning angle (degrees) 

Oblique shock angle (degrees) 

Mach Angle (degrees) 

Prandtl-Meyer angle (degrees) 

Static mass density 
Mass density 


Subscripts 

* 

1 

2 


Description 

Critical flow condition (i.e., conditions where the local fluid 
velocity is equal to the local speed of sound) 

Upstream flow property 
Downstream flow property 
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3. Quick Reference Tables 



PROPERTIES OF ISENTROPIC FLOW, PRANDTL-MEYER FLOW, 
AND NORMAL SHOCKS 

ames 

Solves the equations for isentropic flow, Prandtl-Meyer flow, and 
normal shocks to obtain flow properties. 

amesplt 

Plots the properties for isentropic flow, Prandtl-Meyer flow, and the 
normal shocks as a function of Mach number. 

ameserr 

Consistency check for function ames. Computes and plots, as a function 
of Mach number, errors in ames calculations. 

isentbl 

Generates text file containing a table of the isentropic flow properties. 

nshktbl 

Generates text file containing a table of Prandtl-Meyer flow and normal 
shock properties. 


PROPERTIES OF OBLIQUE SHOCKS 

oblqshck 

Solves the oblique shock equations for both weak and strong shock 
angles. 

obi qwl2 

Solves the oblique shock equations to obtain downstream flow 
properties as a function of upstream flow properties. 

oblqw21 

Solves the oblique shock equations to obtain upstream flow properties 
as a function of downstream flow properties. 

deltason 

Computes the theoretical deflection angle that reduces supersonic flow 
to sonic conditions. 

deltamax 

Computes the theoretical maximum angle through which supersonic 
flow may be deflected or turned without separation. 


PROPERTIES OF FANNO-LINE FLOW 

fanno 

Solves the Fanno line equations to obtain properties of flow with 
friction. 

fannoplt 

Plots the Fanno line flow properties as a function of Mach number. 

fannoerr 

Consistency check for function fanno. Computes and plots, as a 
function of Mach number, errors in fanno calculations. 

fannotbl 

Generates text file containing a table of the Fanno-line flow properties. 
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PROPERTIES OF RAYLEIGH-LINE FLOW 

rayleigh 

Solves the Rayleigh-line equations to obtain properties of flow heating 
or cooling. 

raylplt 

Plots the Rayleigh-line flow properties as a function of Mach number. 

raylerr 

Consistency check for function rayleigh. Computes and plots, as a 
function of Mach number, errors in rayleigh calculations. 

rayltbl 

Generates text file containing a table of the Rayleigh-line flow 
properties. 
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4. Function Reference Guide 


ames 


Purpose 

Solve the equations for isentropic flow, both subsonic and supersonic, Prandtl-Meyer 
expansion, and normal shocks. 

Synopsis 

ames 

Properties=ames( Varln , Val uesln , VarsOut) 

Properties=ames( Varln , Val uesln , VarsOut , Gamma) 

[Properties , PI tLbl s]=ames (Varln , Val uesln , VarsOut , Gamma) 

Description 

ames by itself calls amespl t which plots normalized versions of the isentropic flow, 
Prandtl-Meyer, and normal shock functions versus Mach number. 

Properti es=ames( Varln , Val uesln , VarsOut) , given a number designating one ofthe 
flow properties listed in Table 4.1 and a value or vector of values for that flow property, 
ames computes corresponding values for isentropic flow, Prandtl-Meyer flow, and 
normal shock functions. Varln is a scalar that specifies the property used as the input 
(independent variable). Val uesln may be a scalar or vector and contains values of the 
independent variable for which the other properties will be computed. VarsOut contains a 
list of Indices corresponding to the flow properties listed in Table 4.1. Indices specified 
in VarsOut may be in any order and may be repeated as desired by the user. Results are 
returned in the Properties matrix. Columns in this matrix correspond to indices 
specified in VarsOut. Rows ofthe Properties matrix contain results corresponding to 
the elements of Val uesln. 

Note that, when properties 5, 6, or 7 are used as the independent variable, the solution is 
double- valued. The double-valued solution is provided by making Properties a cell 
array. Properti es{l} contains values of the solution associated with the smaller Mach 
number, while Properti es{2) contains the solution associated with the larger Mach 
number. 

Properti es=ames( Varln , Val uesln , VarsOut , Gamma) provides a mechanism for 
specifying values for the ratio of specific heats of the working fluid via Gamma. Gamma is 
optional. If unspecified, a value of 1.4, the value of the ratio of specific heats of air at 
standard temperature and pressure, is used. If specified. Gamma may be defined as either a 
scalar or a vector. If it is a vector, it must have the same length as Val uesln. 

[Properties , PI tLbl s]=ames (Varln , Val uesln , VarsOut , Gamma) , in addition to 
returning the properties of the fluid at user specified conditions, also returns a cell array, 
PI tLbl S, containing text strings that may be used when plotting the results. 
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Table 4.1 — Description of Flow Properties Computed by Function ames 


REF. 

INDEX 

PROPERTY 

REF. 1 

DESCRIPTION 

ISENTROPIC FLOW PROPERTIES (VALID FOR ALL M): 

1. 

M or Mj 


Mach number 

2. 

P/P, 

Eq. 44 

Ratio of static to total pressure 

3. 

p/p t 

Eq. 45 

Ratio of static to total density 

4. 

m 

Eq. 43 

Ratio of static to total temperature 

5. 

p 

Pg- 1 

iK > 

6. 

q/P, 

Eq. 48 

Ratio of dynamic to total pressure 

7. 

A/ A" 

Eq. 80 

Ratio of flow area to critical flow area 

8. 

v/v * 

Eq. 50 

Ratio of flow velocity to critical flow velocity 

PRANDTL-MEYER FLOW (VALID FOR M>1): 

9. 

V 

Eq. 171 

Prandtl-Meyer angle (degrees) 

10. 

F 

Pg- 1 

Mach Angle ( degrees), sin - 1 (l / M ) 

NORMAL SHOCK PROPERTIES (VALID FOR M>1): 

11. 

m 2 

Eq. 96 

Mach number downstream of a normal shock 

12. 

Pl/Pl 

Eq. 93 

Static pressure ratio across a normal shock 

13. 

P2 /Pi 

Eq. 94 

Static density ratio across a normal shock 

14. 

TJT X 

Eq. 95 

Static temperature ratio across a normal shock 

15. 

Pi/Pa 

Eq. 99 

Total pressure ratio across a normal shock 

16. 

P/P,, 2 

Eq. 100 

Ratio of static pressure upstream of a normal shock 
to total pressure downstream of the same shock 


Algorithm 

ames determines the desired flow properties by first obtaining a Mach number solution 
for each value, Val uesln, of the user specified flow property, Varln. These Mach 
numbers are then used to compute the other properties, VarsOut, specified by the user. 
Most of the flow equations may be manipulated analytically to obtain Mach number as a 
function of the other properties. However, some nonlinear relationships exist which have 
no simple analytical solution. In these cases, MATLAB’s fminbnd function is used 
determine an approximate solution for Mach number from the nonlinear equations. The 
search is arbitrarily constrained to Mach numbers less than 100. Solutions associated with 
Mach numbers larger than 100 are returned as NaN (i.e., not a number). 


See Also 

ameserr, amespl t, 1 sentbl , and nshcktbl 


Example 4.1: 

Determine the properties of air at Mach 2. 
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» ames (1,2,1:16) 
ans = 

Columns 1 through 5 


2.0000 

0.1278 

0 . 

.2300 

0.5556 

1 . 

.7321 

Columns 

6 through 10 






0.3579 

1.6875 

1 . 

.6330 

26.3798 

30. 

.0000 

Columns 

11 through 15 






0.577 

44.5000 

2. 

.6667 

1.6875 

0 . 

.7209 


Column 16 
0.1773 


Example 4.2: 

Given a normal shock with downstream Mach number of 0.85, determine the Mach 
number upstream of the shock. 

» ames(ll, 0.85,1) 
ans = 

1.1876 

Example 4.3: 

Determine the properties of air when A/ A* = 3.007. 

» properties=ames(7,3.007,l:16) 
properties = 

[1x16 double] [1x16 double] 

» properties(l) 
ans = 

Columns 1 through 5 


0.1970 

0.9733 

0.9809 

0.9923 

0.9804 

Columns i 

6 through 10 




0.0264 

3.0070 

0.2149 

NaN 

NaN 

Columns 

11 through 15 




NaN 

NaN 

NaN 

NaN 

NaN 


Column 16 
NaN 

» properties{2) 
ans = 

Columns 1 through 6 


2.6399 

0.0471 

0.1128 

0.4177 

2. 

.4432 

Columns i 

6 through 10 





0.2299 

3.0070 

1.8691 

42.3049 

22. 

.2597 

Columns 

11 through 15 





0.5005 

7.9638 

3.4935 

2.2796 

0 . 

.4453 


Column 16 
0.1058 
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Example 4.4: 

Plot the Mach number downstream of a normal shock as a function of the Mach number 
upstream of the shock. 

Ml=l : 0 . 1 : 10 : 

[M2,Lbls]=ames(l,Ml,ll) ; 
pi ot (Ml . M2) ; 

xl abel ( ' M_l' ) ; ylabel (Lbls{l}) : 



Figure 4.1 . — Result of using function ames to compute 
Mach number variations across a normal shock. 
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ameserr 


Purpose 

Show the computational errors that result when using function ames to solve the 
equations for isentropic flow, Prandtl-Meyer expansion, and normal shocks. 

Synopsis 

ameserr 

[error ,Ml]=ameserr 

Description 

ameserr computes the error between Mach numbers used as inputs to function ames and 
Mach numbers calculated from the output of function ames. The results are plotted as 
absolute and percent errors versus Mach number for each of the flow functions shown in 

Table 4.1. 

[error , Ml]=ameserr returns the computed error in error. If specified. Ml contains the 
initial vector of Mach numbers. 

Algorithm 

ameserr first generates a logarithmically spaced vector of 250 Mach numbers from 0.01 
to 10. This vector also includes critical Mach number values where numerical stability is 
important, such as saddle points, ameserr then uses function ames to calculate each of 
the isentropic flow properties and the normal shock properties corresponding to those 
Mach numbers. The functions of Mach number, obtained from ames, are then used as 
input to the ames function in order to obtain a Mach number which corresponds to the 
function value. Theoretically, the initial and computed Mach numbers should be the 
same. In general, they are not due to round off, truncation, convergence, and/or 
optimization errors. The difference in the two Mach numbers is returned as the error in 
the calculations. 

See Also 

ames, amesplt, isentbl, and nshcktbl 


Example 4.5: 

Compute and plot the errors the errors that result from running ameserr. Plots are shown 
in Figure 4.2(a to g). 

» ameserr 
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1 


Test Consistency of AMES.M 



(b) 


Figure 4.2. — Output of function ameserr as computed on 
an Intel Pentium4 processor-based computer running 
MATLAB® 7. 
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Test Consistency of AMES.M 



x ICf 8 



Mach No. 


(c) 



(d) 

Figure 4.2. — Output of function ameserr as computed on an 
Intel Pentium4 processor-based computer running 
MATLAB® 7 (continued). 
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Test Consistency of AMES.M 



Mach No. 


(e) 


Test Consistency of AMES.M 



A 


i l 

Mach No. 


(f) 


Figure 4.2. — Output of function ameserr as computed on 
an Intel Pentium4 processor-based computer running 
MATLAB® 7 (continued). 
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LU 


Q. 

"'M 


Q. 


O 

LU 


O 

LU 


Mach No. 


(g) 


Q_ 

"Am 

oT 

o 

uj 


Q_ 

oT 


o 

LU 


Test Consistency of AMES.M 



x 10" 7 



Mach No. 



(h) 

Figure 4.2. — Output of function ameserr as computed on an 
Intel Pentium4 processor-based computer running 
MATLAB® 7 (continued). 
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amesplt 


Purpose 

Plots normalized properties for isentropic flow, Prandtl-Meyer expansion, and normal 
shocks as a function of Mach number. 

Synopsis 

amesplt 

amespl t ( MNmi n , MNmax) 
amespl t ( MNmi n , MNmax , Npts) 
amespl t ( MNmi n , MNmax , Npts , Gamma ) 

Description 

amespl t uses function ames to compute and plot the isentropic and normal shock flow 
properties at 250 points between Mach 0.01 and Mach 10 when the ratio of specific heats 
of the fluid is 1.4. 

amespl t( MNmi n , MNmax) plots results for a range of user specified Mach numbers where: 
MNmi n is the minimum Mach number; and MNmax is the maximum Mach number. 

amespl t ( MNmi n , MNmax, Npts) in addition to allowing the user to specify the range of 
Mach numbers used, this form allows the user to specify the number of data points, Npts, 
used to plot each curve. 

amespl t( MNmi n, MNmax, Npts , Gamma) in addition to allowing the user to specify Mach 
number, and number of points per curve, this form also allows the user to specify a scalar 
value for the ratio of specific heats, Gamma, of the fluid. 

Algorithm 

amespl t first generates a logarithmically spaced vector of 250 Mach numbers from 0.01 
to 10. This vector also includes critical Mach number values where numerical stability is 
important, such as solution saddle points, amespl t then uses this vector as inputs to 
function ames which is used to calculate each of the isentropic flow properties and the 
normal shock properties corresponding to those Mach numbers. The resulting values are 
normalized and plotted versus Mach number to provide the user a graphical 
understanding of the relationship between flow properties and Mach number. 

See Also 

ames, amesplt, isentbl, and nshcktbl 


Example 4.6: 

Plot normalized isentropic flow and normal shock properties as a function of Mach 
number. The resulting plots are shown in Figure 4.3 (a and b). 

» amesplt 
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AMES.M: Isentropic Functions, Table Columns 2-8 



(a) 

AMES.M: Normal Shock Functions, Table Columns 9-16 



(b) 


Figure 4.3. — Normalized isentropic and normal shock 
functions as generated by function amesplt. 
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deltamax 


Purpose 

For steady state supersonic flow with compressive turning, del tamax computes the 
maximum flow deflection angle (8) that can occur without producing separation of the 
flow from the turning surface. Also, optionally calculates the angle of the oblique shock 
(0) that results from turning the flow. Both angles have units of degrees. See Figure 4.4 
for a graphical representation of the flow situation. 

Synopsis 

del tamax 

Delta=deltamax(Ml) 

[Del ta,Theta]=deltamax(Ml, Gamma) 

Description 

del tamax by itself, computes and plots the maximum flow deflection and resulting 
oblique shock angle for a range of Mach numbers from 1.0 to 15. 

Del ta=del tamax(Ml) computes and returns the maximum flow deflection angle, del ta, 
in degrees for user specified Mach numbers. Ml. Ml may be a scalar, vector, or matrix. 

[Del ta,Theta]=del tamax(Ml , Gamma) uses optional input Gamma, the ratio of specific 
heats for the working fluid, to calculate the turning angle Del ta and additionally the 
angle, Theta, of the oblique shock that results from turning the flow. Gamma has a 
default value of 1.4 and must be a scalar or have dimensions equivalent to Ml. The 
dimensions of Del ta and Theta, and the values therein, correspond to the dimensions of 
Ml. 

Algorithm 

If no input parameters are specified by the user, del tamax first generates a vector of 
upstream Mach numbers. The function then uses the Mach number(s) to calculate the 
maximum angle, 0 max , of an oblique shock that can occur without separation. The shock 
angle is then used with the Mach number(s) to calculate the associated flow deflection 
angle, &max . 


r 

/ 

/ 

Shock wave — 

/ Region 2 



Figure 4.4. — Oblique shock diagram. 
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The equation used to calculate 0 max is: 


6max — s i n 1 . 


1 Jfr + 1 W 


yMf 



1 + 

2 1 16 1 


-1 


(4.1) 


The equation used to calculate 5 max is: 


8 


max 


= tan -1 


(Mj- sin- 9 max 1 )c ot 0 max 

|(Y + l) M i 2 - M \ sin 2 Omax + 1 


(4.2) 


Similar equations may be found in reference 1, pp. 9 and 12; (ref. 2), p. 586; and (ref. 4), 
pp. 315 to 316. 

See Also 

del tason, obi qshck, obi qw!2, and obi qw21 


Example 4.7: 

Calculate and plot the maximum compressive turning angle and oblique shock angle for 
airflow over a range of Mach numbers from 1 to 15. 

» deltamax 


Maximum Deflection Angle vs. Mach Number 




Figure 4.5. — Results of function deltamax 
showing maximum turning angle and the angle 
of the resulting oblique shock as a function of 
upstream Mach number. 
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Example 4.8: 

Calculate the maximum compressive turning angle and oblique shock angle for steam 
flowing at Mach numbers from 1.5 to 3.0. The ratio of specific heats for steam is 1.327 at 
standard temperature. 

» [Del ta ,Theta]=del tamax( 1 .5:0. 1:3. 0,1. 327) 

Delta = 

Columns 1 through 5 


12.6726 

15.3598 

17.8660 

20. 

.1780 

22. 

.2960 

Columns i 

5 through 10 






24.2282 

25.9869 

27.5862 

29. 

,0404 

30. 

.3634 

Columns : 

11 through 15 





31.5684 

32.6673 

33.6712 

34. 

.5896 

35. 

.4314 


Column 16 
36.2042 

Theta = 

Columns 1 through 5 


66.7820 

66.0774 

65. 

.6264 

65.3536 

65. 

.2072 

Columns i 

6 through 10 






65.1509 

65.1587 

65. 

.2116 

65.2959 

65. 

.4013 

Columns 

11 through 15 





65.5206 

65.6483 

65. 

.7803 

65.9137 

66. 

.0466 


Column 16 
66.1772 
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deltason 


Purpose 

For steady state supersonic flow with compressive turning, del tamax computes the flow 
deflection angle (8) that results in sonic flow downstream of the resulting oblique shock 
(i.e., M 2 = 1). Also, optionally calculates the angle of the oblique shock (0) that produces 
sonic flow. Both angles have units of degrees. See Figure 4.4 (pp. 19) for a graphical 
representation of the flow situation. 

Synopsis 

deltason 

Delta=deltason(Ml) 

[Del ta , Theta]=del tason ( Ml , Gamma ) 

Description 

del tason by itself, computes and plots the sonic flow deflection angle and the resulting 
oblique shock angle for a range of Mach numbers from 1.0 to 15. 

Del ta=del tason (Ml) computes and returns the flow deflection angle, Del ta, that results 
in sonic flow downstream. Values of Del ta are in degrees for user specified Mach 
numbers, Ml. Ml may be a scalar, vector, or matrix. 

[Del ta,Theta]=del tason (Ml , Gamma) uses optional input Gamma, the ratio of specific 
heats for the working fluid, to calculate the turning angle Del ta and additionally the 
angle, Theta, of the oblique shock that results from turning the flow. Gamma has a 
default value of 1.4 and must be a scalar or have dimensions equivalent to Ml. The 
dimensions of Del ta and Theta, and the values therein, correspond to the dimensions of 
Ml. 

Algorithm 

If no input parameters are specified by the user, del tason first generates a vector of 
upstream Mach numbers. The function then uses the Mach number(s) to calculate the 
angle, 0», of an oblique shock that produces sonic flow downstream of the shock. The 
shock angle is then used with the Mach number(s) to calculate the associated flow 
deflection angle, 8*. 

The equation used to calculate 0* is: 


0* = sin 


1 yA/j 2 


(y + l)Mj2 -(3-y) 


' J(' Y + 1 ) 


9+ y Vx M 2 + (y±1) M 4' 

o i i z: 1 


16 


16 


(4.1) 


The equation used to calculate 8* is: 
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(4.2) 


8. = tan- 1 , K 2s '° 29 --» ot9 - 

i(Y + l)M 1 2 -M 1 2 sin 2 0* + l 

Similar equations may be found in reference 1, pp. 9 and 12; (ref. 2), p. 586; and (ref. 4), 
pp. 315 to 316. 

See Also 

del tamax, obi qshck, obi qw!2, and obi qw21 


Example 4.9: 

For airflow over a range of Mach numbers from 1 to 15, calculate and plot the 
compressive turning angle and associated oblique shock angle that results in sonic flow 
downstream of the shock. 

» deltason 


Deflection Angle for Sonic Flow vs. Mach Number 




Figure 4.6. — Results of function deltason 
showing sonic turning angle and the angle of 
the resulting oblique shock as a function of 
upstream Mach number. 


Example 4.10: 

Given a flow of hydrogen gas at a several Mach numbers from 1.0 to 2.0, at each Mach 
number, calculate the compressive turning angle that produces sonic flow and the 
associated oblique shock angle. The ratio of specific heats for hydrogen is 1.667 at 
standard temperature. 

» format short e 

» [Del ta , Theta]=del tason (1.0:0.1:2.0,1.667) 


NASA/TM— 2006-214086 


24 


Delta = 

Columns 1 through 4 
-1.4216e-022 1.2526e+000 3.2647e+000 
Columns 5 through 8 
7.8286e+000 1.0071e+001 1.2192e+001 
Columns 9 through 11 
1.5967e+001 1.7612e+001 1.9103e+001 


Theta = 

Columns 1 through 4 

9.0000e+001 7.3209e+001 6.7949e+001 

Columns 5 through 8 

6.2929e+001 6.1691e+001 6.0910e+001 

Columns 9 through 11 

6.0164e+001 6.0034e+001 5.9998e+001 
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5.5235e+000 

1.4161e+001 

6.4868e+001 

6.0434e+001 




fanno 


Purpose 

Solve the equations for one-dimensional steady adiabatic flow in a constant area duct 
with friction. 

Synopsis 

Properties=fanno( Varln , Val uesln , VarsOut) 

Properties=fanno( Varln , Val uesln , VarsOut , Gamma) 

[Properties , PI tLbl s]=fanno( Varln , Val uesln , VarsOut , Gamma) 
fanno 

Description 

Properties=fanno( Varln, Val uesln, VarsOut) , given a number designating one of 
the flow properties listed in Table 4.2 and a value or vector of values for that flow 
property, fanno computes corresponding values for adiabatic frictional flow. Varln is a 
scalar that specifies the property used as the input (independent variable). Val uesln may 
be a scalar or vector and contains values of the independent variable for which the other 
properties will be computed. VarsOut contains a list of Indices corresponding to the flow 
properties listed in Table 4.2. Indices specified in VarsOut may be in any order and may 
be repeated as desired by the user. Results are returned in the Properti es matrix. 
Columns in this matrix correspond to indices specified in VarsOut. Rows of the 
Properties matrix contain results corresponding to the elements ofValuesIn. 

Note that, when properties 4, 6, or 7 are used as the independent variable, the solution is 
double-valued. The double-valued solution is provided by making Properties a cell 
array. Properti es{l} contains values of the solution associated with the smaller Mach 
number, while Properti es{2) contains the solution associated with the larger Mach 
number. 

Properti es=fanno( Varln , Val uesln , VarsOut , Gamma) provides a mechanism for 
specifying values for the ratio of specific heats of the working fluid via Gamma. Gamma is 
optional. If unspecified, a value of 1.4, the value of the ratio of specific heats of air at 
standard temperature and pressure, is used. If specified. Gamma may be defined as either a 
scalar or a vector. If it is a vector, it must have the same length as Val uesln. 

[Properties , PI tLbl s]=fanno( Varln , Val uesln , VarsOut , Gamma) , in addition to 
returning the properties of the fluid at user specified conditions, also returns a cell array, 
PI tLbl S, containing text strings that may be used when plotting the results. 

fanno by itself calls fannopl t which plots the Fanno-line flow properties versus Mach 
number. 
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Table 4.2. — Description of Flow Properties Computed by Function fanno 


REF. 

INDEX 

PROPERTY 

REF. 4 

DESCRIPTION 

1 . 

M ox Mi 


Mach number 

2. 

t/t , ; 

Eq. 5.31 

Ratio of static temperature at Mi to static temperature 
at sonic conditions. 

3. 

p/p* 

Eq. 5.30 

Ratio of static pressure at M t to static pressure at sonic 
conditions. 

4. 

Pt/Pt * 

Eq. 5.34 

Ratio of total pressure at Mi to total pressure at sonic 
conditions. 

5. 

v/v* 

p*/p 

Eq. 5.29 

Ratio of flow velocity at Mi to flow velocity at sonic 
conditions. Also, ratio of static density at sonic 
conditions to static density at Mi. 

6. 

I/I* 

Eq. 3.42 

Ratio of the impulse function at Mi to the impulse 
function at sonic conditions. 

7. 

4JL*/ D H 

Eq. 5.35 

Friction factor 


Algorithm 

fanno determines the desired flow properties by first obtaining a Mach number solution 
for each value, Val uesln, of the user specified flow property, Varln. The resulting Mach 
numbers are then used to compute the other properties, VarsOut. Some of the flow 
equations may be manipulated analytically to obtain Mach number as a function of the 
other properties. However, some nonlinear relationships exist which have no simple 
analytical solution. In these cases, MATLAB’s fminbnd function is used determine an 
approximate solution for Mach number from the nonlinear equations. The search is 
arbitrarily constrained to Mach numbers less than or equal to 100. Solutions associated 
with Mach numbers larger than 100 are returned as NaN (i.e., not a number). 


See Also 

fannoerr, fannoplt, and fannotbl 


Example 4.11: 

For air flowing at Mach 3.5, determine the Fanno-line flow properties. 

» fannod.3. 5,1:7) 
ans = 

Columns 1 through 5 

3.5000 0.3478 0.1685 6.7896 2.0642 

Columns 6 through 7 
1.2743 0.5864 


Example 4.12: 

For a range of friction factors from 0.5 to 1.0, determine the static pressure ratio (P/P*) 
and upstream Mach number of air flowing adiabatically through a constant area duct. 

» fric=[0. 5:0. 1:1.0] ' ; properties=fanno(7, fric, [3,1]) 
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properties = 

[6x2 double] [6x2 double] 

» [fric properties{l}] 
ans = 


0.5000 

1.7706 

0.5977 

0.6000 

1.8459 

0.5748 

0.7000 

1.9154 

0.5551 

0.8000 

1.9804 

0.5378 

0.9000 

2.0416 

0.5225 

1.0000 

2.0996 

0.5087 


» [fric properties{2}] 
ans = 


0.5000 

0.2359 

2.8603 

0.6000 

0 . 1583 

3.6302 

0.7000 

0.0850 

5.1405 

0.8000 

0.0148 

12.7693 

0.9000 

NaN 

NaN 

1.0000 

NaN 

NaN 


Here, properties(l) is the subsonic solution, and properties{2} is the supersonic 
solution. Also, column 1 of ans contains values for friction factor. Column 2 and 3 
contain corresponding solutions for pressure ratio and Mach number, respectively. Note 
that NaN results imply solution has exceeded internal limit for intermediate Mach number 
calcualtion (i.e., M > 100). 
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fannoerr 


Purpose 

Show the computational errors that result when using function fanno to solve equations 
for one-dimensional steady adiabatic flow in a constant-area duct with friction. 

Synopsis 

fannoerr 

[error , Ml]=fannoerr 

Description 

fannoerr computes the error between Mach numbers used as inputs to function fanno 
and Mach numbers calculated from the output of function fanno. The results are plotted 
as absolute and percent errors versus Mach number for each of the flow functions shown 

in Table 4.2. 

[error , Ml]=fannoerr returns the computed error in error. If specified. Ml contains the 
initial vector of Mach numbers. 

Algorithm 

fannoerr first generates a logarithmically spaced vector of 250 Mach numbers from 0.01 
to 10. This vector also includes critical Mach number values where numerical stability is 
important, such as saddle points, fannoerr then uses function fanno to calculate each of 
the Fanno-line flow properties corresponding to those Mach numbers. The functions of 
Mach number, obtained from fanno, are then used as input to the fanno function in order 
to obtain a Mach number which corresponds to the function value. Theoretically, the 
initial and computed Mach numbers should be the same. In general, they are not due to 
round off, truncation, convergence, and/or optimization errors. The difference in the two 
Mach numbers is returned as the error in the calculations. 

See Also 

fanno, fannoplt, and fannotbl 

Example 4.13: 

Compute and plot the errors the errors that result from running fannoerr. Plots are 
shown in Figure 4.7(a to d) 

» fannoerr 
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(a) 


Test Consistency of FAN NO. M 



x 10" 7 



Mach No. 


(b) 

Figure 4.7. — utput of function fannoerr as computed on an 
Intel Pentium4 processor-based computer running 
MATLAB® 7. 
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Test Consistency of FAN NO. M 



x 10" 7 



10" 2 10"' 10° 10 1 


Mach No. 


(c) 



Mach No. 


(d) 


Figure 4.7. — utput of function fannoerr as computed on 
an Intel Pentium4 processor-based computer running 
MATLAB® 7 (continued). 
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fannoplt 


Purpose 

Plot properties for Fanno-line flow, i.e., one-dimensional steady adiabatic flow in a 
constant-area duct with friction. 

Synopsis 

fannoplt 

f annopl t ( MNmi n , MNmax) 
fannopl t ( MNmi n , MNmax , Npts ) 
fannopl t (MNmi n , MNmax , Npts , Gamma) 

Description 

fannopl t uses function fanno to compute and plot the Fanno-line flow properties at 
250 points between Mach 0.05 and Mach 2.5 when the ratio of specific heats of the fluid 
is 1.4. This plot resembles Figure 5.4 in (ref. 4). 

fannopl t( MNmi n , MNmax) plots results for a range of user specified Mach numbers 
where: MNmi n is the minimum Mach number; and MNmax is the maximum Mach number. 

fannopl t ( MNmi n , MNmax, Npts) in addition to allowing the user to specify the range of 
Mach numbers used, this form allows the user to specify the number of data points, Npts, 
used to plot each curve. 

fannopl t( MNmi n, MNmax, Npts , Gamma) in addition to allowing the user to specify Mach 
No. and number of points per curve, this form also allows the user to specify a scalar 
value for the ratio of specific heats, Gamma, of the fluid. 

Algorithm 

fannopl t first generates a logarithmically spaced vector of 250 Mach numbers from 0.05 
to 2.5. This vector also includes critical Mach number values where numerical stability is 
important, such as solution saddle points, fannopl t then uses this vector as inputs to 
function fanno which is used to calculate each of the isentropic flow properties and the 
normal shock properties corresponding to those Mach numbers. The resulting values are 
plotted versus Mach number to provide the user a graphical understanding of the 
relationship between flow properties and Mach number. 

See Also 

fanno, fannoerr, and fannotbl 

Example 4.14: 

Plot Fanno-line flow properties over a range of Mach numbers from 0.05 to 2.5. The 
resulting plot is shown in Figure 4.8. 

» fannoplt 
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FANNO.M: Fanno Flow Functions 



Figure 4.8. — Fanno-line flow properties as generated by function fannoplt. 
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fannotbl 


Purpose 

Generate a text file containing tables of Fanno-line flow properties. The tables generated 
by this function may be useful when computational solution of the equations is not 
practical. 

Synopsis 

fannotbl 

fannotbl ( Fi 1 ename , Mn , Gamma ) 

Description 

fannotbl uses function fanno to generate a table of values for Fanno-line flow 
properties as a function of Mach numbers from 0.01 to 10. Properties 2 through 7 of 
Table 4.2 are written to the text file, fannotbl . txt. 

fannotbl (Fi 1 ename , Mn , Gamma) computes the flow functions and writes the ASCII data 
to the file specified by the string variable, Fi 1 ename. Functions are evaluated at Mach 
numbers specified in Mn. Gamma is an optional scalar variable specifying the ratio of 
specific heats of the working fluid. If unspecified, a value of 1 .4 is used for Gamma . 

See Also 

fanno, fannoplt, and fannotbl 


Example 4-15: 

Create a table containing values for Fanno-line flow functions over a range of Mach 
numbers from 0.50 to 0.70 in increments of 0.01. Results are shown in Table 4.3 on the 
following page. 

» fannotbl ( ' fannotbl . txt .',0.5:0.01:0.7) 
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Table 4.3. — Output of function fannotbl for a range of Mach numbers from 0.5 to 0.7 
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isentbl 


Purpose 

Generate a text file containing tables of isentropic flow properties. The tables generated 
by this function may be useful when computational solution of the equations is not 
practical. 

Synopsis 

isentbl 

i sentbl ( Fi 1 ename , Mn , Gamma ) 

Description 

isentbl uses function ames to generate a table of values for isentropic flow properties as 
a function of Mach numbers from 0.01 to 10. Properties 2 through 8 of Table 4.1 are 
written to the text file, i sentbl . txt. 

i sentbl ( Fi 1 ename , Mn , Gamma ) computes the flow functions and writes the ASCII data 
to the file specified by the string variable, Fi 1 ename. Functions are evaluated at Mach 
numbers specified in Mn. Gamma is an optional scalar variable specifying the ratio of 
specific heats of the working fluid. If unspecified, a value of 1 .4 is used for Gamma . 

See Also 

ames, fannotbl , nshcktbl , and rayl tbl 

Example 4-16: 

Create a table containing isentropic functions for a range of Mach numbers from 0.9 to 
1.1 in increments of 0.01. Results are shown in Table 4.4 on the following page. 

» isentbl ( 'isentbl .txt. ' ,0.9:0.01:1.1) ; 
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Table 4.4. — Output of function isentbl for a range of Mach numbers from 0.9 to 1.1 
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nshktbl 


Purpose 

Generate a text file containing tables of supersonic flow and normal shock properties. 
The tables generated by this function may be useful when computational solution of the 
equations is not practical. 

Synopsis 

nshktbl 

nshktbl (Filename.Mn, Gamma) 

Description 

nshktbl uses function ames to generate a table of values for supersonic flow and normal 
shock properties as a function of Mach numbers from 1 to 10. Properties 9 through 16 of 
Table 4.1 are written to the text file, nshktbl . txt. 

nshktbl (Filename.Mn, Gamma) computes the flow functions and writes the ASCII data 
to the file specified by the string variable, Fi 1 ename. Functions are evaluated at Mach 
numbers specified in Mn. Gamma is an optional scalar variable specifying the ratio of 
specific heats of the working fluid. If unspecified, a value of 1 .4 is used for Gamma . 

See Also 

ames, isentbl, fannotbl, and rayltbl 


Example 4-17: 

Create a table containing supersonic flow and normal shock functions for a range of 
Mach numbers from 1.0 to 2.5 in increments of 0.1. Results are shown in Table 4.5 on 
the following page. 

» nshktbl ( ' nshktbl . txt .',1.0:0. 1:2. 5): 
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Table 4.5. — Output of function nshktbl for a range of Mach numbers from 1.0 to 2.5 
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oblqshck 


Purpose 

For steady state supersonic flow, obi qshck computes the angle of the oblique shock that 
results from compressively turning the flow through angle (8). One of two solutions will 
occur for each Mach number specified, a weak shock solution, or a strong shock solution. 
Angles have units of degrees. The flow situation is similar to that depicted in Figure 4.4 
(PP- 19). 

Synopsis 

oblqshck 

ThetaW=obl qshck(Ml , del ta) 

[ThetaW,ThetaS]=obl qshck(Ml , del ta , gamma) 

Description 

obi qshck by itself generates a plot showing shock angle vs. deflection angle for lines of 
constant Mach from 1 to 20. This plot is a representation of Chart 2 in (ref. 1). 

[ThetaW,ThetaS]=oblqshck(Ml .delta) computes both the weak oblique shock angle, 
ThetaW, and the strong oblique shock angle, ThetaS, that are the result of compressively 
turning supersonic flow at Mach number. Ml, through angle del ta. Ml and del ta may be 
either scalars or vectors. If both are vectors, they must have identical dimensions. 

[ThetaW, ThetaS , DELmax, DELson]=obl qshck(Ml , del ta , gamma) uses optional scalar 
input Gamma to specify the ratio of specific heats of the working fluid. If unspecified, a 
value of 1.4 is used for Gamma. In addition to returning the oblique shock angles, this 
form also returns the maximum flow deflection angle, DELmax, and the flow deflection 
angle that results in sonic flow downstream of the oblique shock, DELson. 

Algorithm 

obi qshck uses the upstream Mach number, the flow deflection angle, and the ratio of 
specific heats to calculate the solution of the cubic equation for both weak and strong 
oblique shock angles using the method given in (ref. 5). If flow deflection angles are 
specified as outputs, obi qshck uses functions del tamax and del tason to compute 
values for those parameters at Mach number(s). Ml. 

See Also 

del tason, del tamax, obi qwl2, and obi qw21 


Example 4-18: 

Replicate the results in Chart 2 of (ref. 1). 

» oblqshck 
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OBLQSHCK: NACA Report 1135, Chart 2 



Figure 4.9. — Oblique shock angle versus deflection angle 
for lines of constant Mach number. 


Example 4-19: 

Given freestream airflow at Mach 2.2, calculate the shock angle that results from turning 
the flow through a range of deflection angles from zero to the maximum deflection 
possible without separating the flow. 

» Ml=2.2, DELmax=deltamax(M) 

Ml = 

2.2000 

DELmax = 

26.1028 

» delta=DELmax*[0:10] 710; Ml=Ml*ones(size(delta) ) ; 

» ThetaW=oblqshck(Ml, delta) ; [delta ThetaW] 
ans = 


0 

27.0357 

2.6103 

29.0843 

5.2206 

31.2899 

7.8308 

33.6666 

10.4411 

36.2338 

13.0514 

39.0200 

15.6617 

42.0707 

18.2719 

45.4660 

20.8822 

49.3676 

23.4925 

54.2056 

26.1028 

64.6203 
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oblqwl2 


Purpose 

For steady state supersonic flow, obi qwl2 uses the upstream fluid properties to compute 
properties of the flow downstream of a weak oblique shock. The flow situation is similar 
to that depicted in Figure 4.4 (pp. 19). 

Synopsis 

obi qwl2 

M2=oblqwl2 (Ml .delta) 

[M2 , Theta , PT rati o]=obl qwl2 ( Ml , del ta , gamma ) 

Description 

obi qwl2 by itself generates a series of plots showing properties of oblique shocks for 
lines of constant Mach from 1 to 20. Figure 4.9(a and b) replicates the weak shock 
portions of Charts 2 and 4 in (ref. 1). Figure 4.9(c) shows variations in total pressure 
across and oblique shock as a function shock angle. In these figures, the solid lines 
represent lines of constant Mach number. 

M2=obl qwl2 (Ml , del ta) computes the Mach number, M2, downstream of an oblique 
shock that results from compressively turning supersonic flow at Mach number, Ml, 
through angle del ta. Ml and del ta may be either scalars or vectors. If both are vectors, 
they must have identical dimensions. 

[M2, Theta, PTratio]=oblqwl2(Ml, delta, gamma) uses optional scalar input Gamma to 
specify the ratio of specific heats of the working fluid. If unspecified, a value of 1.4 is 
used. In addition to returning the downstream Mach number, this form also returns the 
resulting oblique shock angle as Theta, and the total pressure ratio across the shock, 
Pta/Pt, i, as PTratio. 

Algorithm 

Obi qwl2 uses the upstream Mach number, the flow deflection angle, and the ratio of 
specific heats to calculate the solution of the cubic equation for both weak and strong 
oblique shock angles using the method given in reference 5. It then uses equation 131 and 
142 from reference 1 — alternatively, equation 7.31 and 7.25 from reference 4 — to 
calculate the downstream Mach number and total pressure ratio across the oblique shock. 

See Also 

del tason, del tamax, obi qshck, and obi qw21 


Example 4-20: 

Replicate the results in Chart 2 and 4 of reference 1 . Results are shown in Figure 4. 1 0. 

» obl qshck 
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Oblique Shock Angle (deg) 


OBLQW12: Weak Shock Wave Angle vs. Deflection Angle 



(a) 


OBLQW12: Mach Number Parameter vs. Deflection Angle 



(b) 


Figure 4.10. — Plots generated by function oblqwl2. 
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Example 4-21: 

Given flow deflection of freestream airflow at Mach 2.2, calculate the shock angle that 
results from turning the flow through a range of deflection angles from zero to the 
maximum deflection possible without separating the flow. 

» Ml=2.2, DELmax=deltamax(M) 

Ml = 

2.2000 
DELmax = 

26.1028 

» delta=DELmax*[0:10] 710; Ml=Ml*ones(size(delta)) ; 

» ThetaW=oblqshck(Ml, delta) ; [delta ThetaW] 
ans = 

0 27.0357 

2.6103 29.0843 

5.2206 31.2899 

7.8308 33.6666 

10.4411 36.2338 

13.0514 39.0200 

15.6617 42.0707 

18.2719 45.4660 

20.8822 49.3676 

23.4925 54.2056 

26.1028 64.6203 
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oblqw21 


Purpose 

For steady state supersonic flow, obi qw21 uses the downstream fluid properties to 
compute properties of the flow upstream of a weak oblique shock. The flow situation is 
similar to that depicted in Figure 4.4 (pp. 19). 

Synopsis 

oblqw21 

Ml=oblqw21 (M2 .delta) 

[Ml , Theta , PT rati o]=obl qw21 (M2 , del ta , gamma) 

Description 

obi qw21 generates a consistency check for the oblique shock equations used by functions 
obi qshck, obi qwl2, and obi qw21. The absolute error is plotted as a function of flow 
deflection angle for lines of constant Mach number. 

Ml=obl qw21 (M2 , del ta) computes the Mach number, Ml, upstream of an oblique shock 
that is required to achieve downstream Mach number, M2, after compressively turning the 
flow through angle del ta. M2 and del ta may be either scalars or vectors. If both are 
vectors, they must have identical dimensions. 


[Ml, Theta, PTratio]=oblqw21(M2, delta, gamma) uses optional scalar input Gamma to 
specify the ratio of specific heats of the working fluid. If unspecified, a value of 1 .4 is 
used. In addition to returning the downstream Mach number, this form also returns the 
resulting oblique shock angle as Theta, and the total pressure ratio across the shock, 

P u 2 /P t ,i, as PTratio. 

Algorithm 

obi qw21 uses MATLAB’s fmi nbnd function to find a solution to the nonlinear equation 
for downstream Mach number as a function of upstream Mach number, the flow 
deflection angle, and the ratio of specific heats, fmi nbnd solves the nonlinear equation by 
employing an inline function to compute the error between the user specified Mach 
number, M 2 , and a downstream Mach number, M 2 , G ues s 5 that is calculated using function 
obl qwl2 and a guess for the upstream Mach number, Mi, Guess . The search is arbitrarily 
constrained to Mach numbers less than or equal to 100. Solutions associated with Mach 
numbers larger than 100 are returned as NaN (i.e., not a number). 


See Also 

del tason, del tamax, obl qshck, and obl qwl2 


Example 4-22: 

Generate the consistency check for the oblique shock functions. 


» oblw21 
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0BLQW21 : Consistency Check showing Error 
for Lines Constant Upstream Mach 



Figure 4.1 1 . — Consistency check for oblique shock 
functions oblqshck, oblqw!2, and oblqw21. 


Example 4-23: 

Given the values for downstream Mach number, M2, shown below, calculate the upstream 
Mach number required to produce M2 when the flow deflection angle is 0.5 degrees. Also 
calculate the oblique shock angle and the total pressure ratio across the oblique shock. 

» M2=[ 0.9000 1.0000 1.9819 2.9742 3.9624 ... 

61.9798 69.5625 76.8279 83.7768 90.4113]; 

» delta=0.5; 

» [Ml, Theta, PTratio]=oblqw21(M2, delta) ; 

Warning: FANNO: Some Mach number solutions exceed Internal limit (M 

> 100). Solutions set to NaN 

> In oblqw21 at 119 
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» [M2 ( : ) Ml ( : ) Theta(:) PTratio(:)] 
ans = 


0.9000 

NaN 

NaN 

NaN 

1.0000 

1.0490 

77.7993 

1.0000 

1.9819 

2.0000 

30.4029 

1.0000 

2.9742 

3.0000 

19.8116 

1.0000 

3.9624 

4.0000 

14.8009 

1.0000 

61.9798 

70.0000 

1.1719 

0.9500 

69.5625 

80.0001 

1.0766 

0.9288 

76.8279 

90.0000 

1.0038 

0.9037 

83.7768 

99.9999 

0.9468 

0.8750 

90.4113 

NaN 

NaN 

NaN 


Note that a weak oblique shock solution does not exist for the case where M2 = 0.9. 
Therefore, the associated solutions are defined as NaN. Also, because the search space 
limited to Mach numbers less than or equal to 100, a solution for M2 = 90.41 13, Ml = 
1 10, is not computed by obi qw21. 
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rayleigh 


Purpose 

Solve the equations for one-dimensional steady flow in a constant area duct with heat 
transfer. 

Synopsis 

Properties=rayl eigh( Varln , Val uesln , VarsOut) 

Properties=rayl eigh( Varln , Val uesln , VarsOut , Gamma) 

[Properties , PI tLbl s]=rayl eigh( Varln , Val uesln , VarsOut , Gamma) 
rayleigh 

Description 

Properties=rayleigh( Varln, Val uesln, VarsOut) , given a number designating one 
of the flow properties listed in Table 4.6 and a value or vector of values for that flow 
property, rayl eigh computes corresponding values for flow though a constant area duct 
with heat transfer. Varln is a scalar that specifies the property used as the input 
(independent variable). Val uesln may be a scalar or vector and contains values of the 
independent variable for which the other properties will be computed. VarsOut contains a 
list of indices corresponding to the flow properties listed in Table 4.6. Indices specified 
in VarsOut may be in any order and may be repeated as desired by the user. Results are 
returned in the Properties matrix. Columns in this matrix correspond to indices 
specified in VarsOut. Rows of the Properties matrix contain results corresponding to 
the elements of Val uesln. 

Note that, when properties 2, 3, or 5 are used as the independent variable, the solution is 
double- valued. The double-valued solution is provided by making Properties a cell 
array. Properties(l) contains values of the solution associated with the smaller Mach 
number, while Properties{2) contains the solution associated with the larger Mach 
number. 

Properties=rayl eigh( Varln , Val uesln , VarsOut , Gamma) provides a mechanism for 
specifying values for the ratio of specific heats of the working fluid via Gamma. Gamma is 
optional. If unspecified, a value of 1.4, the value of the ratio of specific heats of air at 
standard temperature and pressure, is used. If specified. Gamma may be defined as either a 
scalar or a vector. If it is a vector, it must have the same length as Val uesln. 

[Properties , PI tLbl s]=rayl eigh( Varln , Val uesln , VarsOut , Gamma) , in addition to 
returning the properties of the fluid at user specified conditions, also returns a cell array, 
PI tLbl S, containing text strings that may be used when plotting the results. 

rayl ei gh by itself calls rayl pi t which plots the Rayleigh-line flow properties versus 
Mach number. The resulting plot is similar to figure 6.5 in reference 4. 
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Table 4.6. — Description of Flow Properties Computed by Function rayleigh 


REF. 

INDEX 

PROPERTY 

REF. 4* 
EQN. 

DESCRIPTION 

1 . 

MorM/ 


Mach number 

2. 

TJT 

Eq. 6.22 

Ratio of total temperature at M ) to total temperature at 
sonic conditions. 

3. 

T /T, 

Eq. 6.19 

Ratio of static temperature at Mj to static temperature 
at sonic conditions. 

4. 

P/P, 

Eq. 6.20 

Ratio of static pressure at Mj to static pressure at sonic 
conditions. 

5. 

P./P,,, 

Eq. 6.23 

Ratio of total pressure at Mj to total pressure at sonic 
conditions. 

6. 

v/v. 

Eq. 6.21 

Ratio of flow velocity at Mj to flow velocity at sonic 
conditions. 


Similar equations may also be found in reference 2 on p. 196. 


Algorithm 

rayl ei gh determines the desired flow properties by first obtaining a Mach number 
solution for each value, Val uesln, of the user specified flow property, Varln. The 
resulting Mach numbers are then used to compute the other properties, VarsOut. Most of 
the flow equations may be manipulated analytically to obtain Mach number as a function 
of the other properties. However, in the case of total pressure ratio, a nonlinear 
relationships exists which has no simple analytical solution. In this case, MATLAB’s 
fmi nbnd function is used determine an approximate solution for Mach number from the 
nonlinear equation. The search is arbitrarily constrained to Mach numbers less than or 
equal to 100. Solutions associated with Mach numbers larger than 100 are returned as 
NaN (i.e., not a number). 

See Also 

rayl err, rayl pi t, and rayl tbl 

Example 4.24: 

For air flowing at Mach 0.72 and 2.85, determine the Rayleigh-line flow properties. 

» rayl eigh(l. [0.72 2. 85], 2:6) 
ans = 

0.9221 1.0026 1.3907 1.0376 0.7209 

0.6685 0.3057 0.1940 3.0014 1.5757 

Example 4.25: 

Given sonic conditions at a point in a constant area duct with total temperature of 1000 K 
and total pressure of 300 kPa, find the Mach number and total pressure as the total 
temperature decreases to 800, 600, and 400 K. 

First, divide duct total temperature by sonic total temperature to obtain the ratio of total 
temperatures. 
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» TTratio = [900 750 600]/1000 
TTratio = 

0.9000 0.7500 0.6000 


Next, use function rayl ei gh to compute the Mach number (ref. index 1) and total 
pressure ratio (ref. index 5) for the speficied total temperature ratios. 

» Properties=rayleigh( 2, TTratio, [1 5]) 

Properties = 

[3x2 double] [3x2 double] 

The solution is double-valued. The subsonic solution is returned as Properties{l}, and 
the supersonic solution is returned as Properties{2}. Mach number values are returned 
in column one of each of the solutions. The subsonic and supersonic solutions for Mach 
number are shown below in columns one and two, respectively. Here, the rows 
correspond to the elements of TT rati o, with the first row corresponding to the element 
one. 

» Ml=[Properties{l}(: ,1) Properties{2}( : , 1)] 

Ml = 

0.6884 1.5368 

0.5423 2.2361 

0.4415 3.7749 

Values for total pressure ratio are returned in column two of each of the cells in 
Properties. The subsonic and supersonic solutions are shown below in columns one and 
two, respectively. Here, the rows correspond to the elements of TTratio, with the first 
row corresponding to the element one. The total pressure ratios are multiplied by the 
pressure ratio at sonic conditions (i.e., 300 kPa) to obtain total pressure at the specified 
total temperature conditions. 

» PT=[Properties{2}( : ,2) ; Properties{2}( : ,2)]*300 
PT = 

1.0e+003 * 

0.3139 0.3421 

0.3291 0.5379 

0.3416 2.0328 

Note that NaN results would imply that the intermediate solution for Mach number 
exceeded the internal limit (i.e., M > 100). 
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raylerr 


Purpose 

Show the computational errors that result when using function rayl eigh to solve 
equations for one-dimensional steady adiabatic flow in a constant-area duct with friction. 

Synopsis 

raylerr 

[error ,Ml]=rayl err 

Description 

rayl err computes the error between Mach numbers used as inputs to function rayl eigh 
and Mach numbers calculated from the output of function rayl eigh. The results are 
plotted as absolute and percent errors versus Mach number for each of the flow functions 
shown in Table 4.6. 

[error , Ml]=rayl err returns the computed error in error. If specified. Ml contains the 
initial vector of Mach numbers. 

Algorithm 

rayl err first generates a logarithmically spaced vector of 250 Mach numbers from 0.01 
to 10. This vector also includes critical Mach number values where numerical stability is 
important, such as saddle points, rayl err then uses function rayl eigh to calculate each 
of the Rayleigh- line flow properties corresponding to those Mach numbers. The functions 
of Mach number, obtained from rayl eigh, are then used as input to the rayl eigh 
function in order to obtain a Mach number which corresponds to the function value. 
Theoretically, the initial and computed Mach numbers should be the same. In general, 
they are not due to round off, truncation, convergence, and/or optimization errors. The 
difference in the two Mach numbers is returned as the error in the calculations. 

See Also 

rayl ei gh, rayl pi t, and rayl tbl 


Example 4.26: 

Perform a consistency check on the calculations in function rayl eigh. Plots are shown in 
Figure 4.12(a to c). 

» raylerr 
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%Error (T7T ,) 


1 


Test Consistency of RAYLEIGH. M 


o 

m 


0 


1 

0 

-1 



(a) 



(b) 


Figure 4.12. — Output of function raylerr as computed on 
an Intel Pentium4 processor-based computer running 
MATLAB® 7. 
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(C) 


Figure 4. 1 2. — Output of function rayl err as computed on 
an Intel Pentium4 processor-based computer running 
MATLAB® 7 (continued). 
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raylplt 


Purpose 

Plot properties for Rayleigh-line flow, i.e., one-dimensional steady adiabatic flow in a 
constant-area duct with friction. 

Synopsis 

raylplt 

rayl pi t ( MNmi n , MNmax) 

rayl pi t ( MNmi n , MNmax , Npts) 

rayl pi t ( MNmi n , MNmax , Npts , Gamma ) 

Description 

rayl pi t uses function rayl eigh to compute and plot the Rayleigh- line flow properties 
at 250 points between Mach 0.05 and Mach 2.5 when the ratio of specific heats of the 
fluid is 1.4. This plot resembles figure 6.5 in reference 4. 

rayl pi t ( MNmi n , MNmax) plots results for a range of user specified Mach numbers where: 
MNmi n is the minimum Mach number; and MNmax is the maximum Mach number. 

rayl pi t ( MNmi n , MNmax, Npts) in addition to allowing the user to specify the range of 
Mach numbers used, this form allows the user to specify the number of data points, Npts, 
used to plot each curve. 

rayl pi t ( MNmi n , MNmax, Npts , Gamma) in addition to allowing the user to specify Mach 
No. and number of points per curve, this form also allows the user to specify a scalar 
value for the ratio of specific heats, Gamma, of the fluid. 

Algorithm 

rayl pi t first generates a logarithmically spaced vector of 250 Mach numbers from 0.05 
to 2.5. This vector also includes critical Mach number values where numerical stability is 
important, such as solution saddle points, rayl pi t then uses this vector as inputs to 
function rayl eigh which is used to calculate each of the isentropic flow properties and 
the normal shock properties corresponding to those Mach numbers. The resulting values 
are plotted versus Mach number to provide the user a graphical understanding of the 
relationship between flow properties and Mach number. 

See Also 

rayl eigh, rayl err, and rayltbl 

Example 4.27: 

Plot Rayleigh-line flow properties over a range of Mach numbers from 0.5 to 2.5. The 
resulting plot is shown in Figure 4.13. 

» raylplt 
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RAYLPLT.M: Rayleigh-line Flow Functions 



Figure 4.13. — Rayleigh-line flow properties as generated 
by function raylplt. 
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rayltbl 


Purpose 

Generate a text file containing tables of Rayleigh-line flow properties. The tables 
generated by this function may be useful when computational solution of the equations is 
not practical. 

Synopsis 

rayltbl 

rayl tbl ( Fi 1 ename , Mn , Gamma ) 

Description 

rayl tbl uses function rayl eigh to generate a table of values for Rayleigh- line flow 
properties as a function of Mach numbers from 0.01 to 10. Properties 2 through 6 of 
Table 4.2 are written to the text file, rayl tbl . txt. 

rayl tbl ( Fi 1 ename , Mn , Gamma ) computes the flow functions and writes the ASCII data 
to the file specified by the string variable, Fi 1 ename. Functions are evaluated at Mach 
numbers specified in Mn. Gamma is an optional scalar variable specifying the ratio of 
specific heats of the working fluid. If unspecified, a value of 1 .4 is used for Gamma . 

See Also 

rayl ei gh, rayl pi t, and rayl tbl 


Example 4-28: 

Create a table containing values for Rayleigh-line flow functions over a range of Mach 
numbers from 0.55 to 059 in increments of 0.01, and from 2.0 to 2.4 in increments of 0.1. 
Results are shown in Table 4.7 on the following page. 

» M1=[0 .55: 0 .01:0.59 2. 0:0. 1:2. 4] ; 

» rayl tbl ( ' rayl tbl . txt . ' , Ml ) 
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Table 4.7. — Output of function rayltbl for a range of Mach numbers specified in Example 4.28 
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